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ABSTRACT 

A least squares formulation of the system divu = p, curlu = £ is 
surveyed from the viewpoint of both finite element and finite difference 
methods* Closely related arguments are shown to establish convergence 
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Introduction 


This paper concerns the three dimensional Cauchy-Riemann type equations 


( 1 ) 


divu = p, curLu = in D 


_u*n = 0 on r ; 


D is a bounded domain in T? with boundary T on which ji is the outward 
normal. The functions p are prescribed and satisfy the compatability 
conditons 


(2) / pdTr = 0, div£ =0 in D; 

D 

these express necessary conditons that the overdetermined first-order system 
(1) has a solution. 

The numerical solution of these equations will be studied from both the 
finite element and finite difference points of view. Indeed, the major goal 
of this paper is to show how both approaches rest on very similar 
foundations. In so doing we hope our study may provide a point of contact 
between those familiar with the, largely separate, literature about each 
method • 

In the case of the finite element method convergence estimates will be 
shown to result quite directly from proof techniques already common to the 
finite element literature. In contrast, many of these same techniques have 
played little or no role in the analysis of finite difference schemes and one 
of our principal objectives lies in clarifying their relevence to finite 


difference methods 
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The Cauchy -Riemann type equations (1) were chosen for study because they 
are representative of a type of first-order systems that arise in problems 
from electromagnetics and fluid dynamics. For such systems, arbitrary finite 
differencce and finite element approximations to (1) are generally 
unsuitable. This is certainly true of all but a few finite element schemes 
based on Galerkin formulations [1], while simple finite difference 
approximations can present, among other difficulties, special problems in 
incorporating boundary conditions accurately. 

In Part A notations used in this paper are introduced and the basic 
intergral identity 

(3) /[ |curlu| 2 +|divu| 2 ] = /|gradu| 2 

D D 

is derived. This identity has found use in many applications (see, for 
example, [3] where it is used in a discussion of the Navier-Stokes 
equations). The derivation given in Section A.l differs from (3) in 
appearence in order to highlight the structural properties of (1) when viewed 
as a first-order system. 

The fact that (1) is well posed is an immediate consequence of (3) and 
the Lax-Milgram Theorem. This is described in Section A. 2 both for 
completeness and because it shows the fundamental role the least squares ideas 
can play in discussing overdetermined systems like (1). 

The derivation of a stable and optimally convergent finite element scheme 
is almost immediate from the material developed in Part A. The arguments, 
while standard, are included in Part B for completeness. It is interesting to 
note that, unlike other least squares approximations to first-order systems, 
these do not impose any restrictions on the finite element spaces [2]. 
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In Part C an analogous development is given for the Keller box-scheme 
based, instead, on a least squares summation formulation. Here a summation by 
parts formula is used to obtain a summation identity analogous to (3) and 
results in a proof that the scheme is second order accurate. Of special 
interest is the fact that this development, also, is not restricted to uniform 
grids nor to grids obtained by images of uniform grids under a global mapping 
function. Thus, like the finite element schemes, it can also be used on 
irregular grids subject only to standard geometric constraints. Finally, we 
remark that the difference scheme employed here (c. f. [10] [12]) involves the 
variables on the same, rather than on staggered, grids in contrast to [5], 
[9]. 


Part A 

A Least Squares Formulation 


Al. An Integral Identity 

In the following, = (x^x^jX^) is a point in 3^ = , 

and (ijk) indicates that the indices i,j,k are restricted to an even 
permutation of (123). Thus 

3 

divu = l 3,u , 

i=l 

cur lvi = (9^ - 3 fc u . ) , (ijk). 


gradu = ( 9 i u j) 
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With 


we define 


3 

(u,v) = [gradu:gradv] = l gradu gradv , 

i=l 


and 


II IX II ^ = /(_U,tl)dTT, 
D 


II uH Q = /|u| 2 dTT 
D 


III I ^ = Hull 2 + Hull 2 . 


The system (1) may be written in matrix form as 


(A. 1.1) 


3 

Lu = I A 3 u = _f 
i=l 


where f 


T 

(P,C) 


and 


A 

1 



0 

0 

0 

-1 





Set 

(A. 1.2) 


B ± = A^A , (ijk) . 


It is easy to verify that 


(A. 1.3) 


T 

A i A i = 


i. 



i = 1,2,3, 


where I is the 


3 * 3 identity matrix 
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Recalling the definition (_u,v) = [gradu:grad_v] , integration by parts 
leads to the identity 

(A. 1.4) (Lti)'^Lv = Ou»v) + divcj^ - fi(u)v 

where, if £ = (q^q^q^, 

(A. 1.5) q ± = (ijk) , 

and fl(u) is a vector with components 9 ^ given by 

(A. 1.6) n ± (u) = ( B k a i 3 iH. T+ Vi 9 # T ^ (ijk) * 

T 

Suppose u_ is smooth; since B- = -B , then 9 fu) = 0* Also, in terms 

K. K. 

of the components of u_ and v_ the component of is easily verified 

to have the form 


Let 

(A. 1.7) 
so that 


q i = - (Vj+VkH’ (ijk) 


qi = v i (a j u j +a k u k ) + ujajvj+a,^) 


Qi = q± - (vjU i )+a k (v k u i )), 


Since the expression in brackets is a surface divergence then 
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Integrating (A. 1.4) and employing Gauss' theorem 


£ (Idl) dn + j^a’iid o 

Suppose _u»_n = v*ii = 0 on T ; the preceding remark implies that 

A 

.g/U = 0 on r so that 

(A. 1.9) /(Lii) T Lvd7r = J(u.,v)dTr. 

D D 


A2. Norm Estimates and Uniqueness 

The goal here is to use tha basic identity (A. 1.9) and the Lax-Milgram 
Theorem to show existence and uniqueness for the system 

(A.2.1) Lii = _f in D, u*n = 0 on T 

where L is defined in (A. 1.1). To do this we must first formulate (A.2.1) 
in terms of a bilinear form a(*,*) on an apropriate space V. In particular 
we put 

T 

(A. 2. 2) a(v,w) = / Lv Lwd'ir = /{divv*divw + curlv # curlw}dTT 

D D 

Moreover, we let 

(A. 2. 3) V = {veL (D):v # ii = 0 on T, gradv e L (D)} 

with 

II vll = {/ (v,v)dir} ^ = {/ [gradv: gradv] dTr} ^ 

D D 


(A. 2. 4) 
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This is a Hilbert space with the associated inner product 

(A. 2. 5) <u,v> = / (u_,v_)d'TT = /(gradu.:gradv)d7r . 

D D 

That the bilinear form a(*, # ) satisfies the conditons of the Lax- 
Milgram Theorem [13] follows immediately from the integral identity (A* 1.9), 
which can be written 

(A. 2. 6) ~ /(u>v)dTT. 

D 

Indeed, putting v_ = u_ we obtain 

(A. 2* 7) “ HuJI^, 

while 

(A. 2. 8) | a. Cu. » v,) | < II ii II B vH 

is also clear. 

It follows that given any bounded linear functional G(») on V there 
is a unique _ueV for which 

(A. 2. 9) a(u_,v) = G(v) all veV; 

moreover, 

(A. 2. 10) Hull < UGH . 

Our final task is to choose G(*) so that (A# 2. 9) is equivalent to (A. 2.1) - 

2 ->-2 

(A. 2. 2). Indeed, given peL (D), (D) we put 

T 

G(v) “ /f. Ludu = /{pdivv + curly} dff . 

D D 


(A. 2. 11) 
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Observe that 

(A. 2. 12) |G(v) | < Hf. II 0 ILv 1 0 = Ilf. II Q II v II 

so that 

(A.2.13) 161 < llf.» 0 = {BpIIq+II^IIq} 1/2 . 

Moreover, (A. 2. 9) is equivalent to 

(A. 2. 14) / |Lu.-f I 2 dTT = /{ |divu-p | 2 +| curlu-t I 2 }d7T 

D D 

be minmized over veV. Thus, if the data p,£ satisfy the compatability 
conditons 


(A.2.15) 


/pdTT = 0, div£ = 0, 
D 


then the min in (A. 2. 14) will be zero and the minimizing function u_eV will 
satisfy (A. 2.1). 

In conclusion, it follows that if 


(A. 2. 16) 


peL^(D) = 


{ ipeL 2 (D) : J^dTr=o} 
D 


(A. 2. 17) 




= {veL 2 (D): divv=0eD} 


are given then there is a unique _ueV such that (A. 2.1) holds. Moreover, 


Hull < ( IIpHq+H ^" 2 ) 1 ^ 2 . 


(A. 2. 18) 
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Part B 

A Finite Element Treatment 


Bl. Least Squares Formulation 

Since the infinite dimensional problem (A. 2*1) - (A. 2*2) has a natural 
characterization (via the Lax-Milgram Theorem) in terms of a least squares 
formulation, it is reasonable to consider approximations based on these 
ideas. Indeed, Let 

(B.l.l) V h “ V 

be a finite dimensional subspace parameterized by h > o. We seek a u^eV^ 
which minimizes 

(B. 1.2) /|Lv h -f| 2 = /{|divv h -p| 2 +|curlv h -£| 2 }cnr 

D D 

as v h varies over V^. Observe that if a(*,*) and G( # ) are defined as 
in Section A. 2, then u^ is a minimizing function if and only if 

(B.1.3) = all v^ey^. 

Moreover, application of the Lax-Milgram Theorem to shows u^ey^ is 

uniquely determined by (B.1.3). Once a basis is chosen for V^, (B.1.3) 

reduces to a set of symmetric positive definite algebraic equations. 
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B2. Error Estimates 

Combining (A. 2. 9) and (B.1.3) we see that 
(B.2.1) a(ta-u_^,v^ ) = 0 all v^ev « 

This orthogonality condition is central to all error estimates. Indeed, first 
note that if u£V^ is given, then (B.2.1) gives 

(B.2.2) a ( u-u h , u-u h ) = a(u_-u^,u-u^) + a(u^-u^,u^-u^) . 

Thus 

~h 2 2 ~h 2 

lu-u 1 = llu-ujl + lllL-U H . 

h — h — 

It follows that is a best approximation in the sense that 

(B.2.3) Hu-u^H = inf llu-u h H < inf flu-u^, 

where the inf is taken over all u 1 in V, • 

— n 

In particular, if consist of piecewise linear elements, then (B.2.3) 

gives 

(B.2.4) Hu-uJI < ChllulL, 

where h is a generic mesh spacing. Here II • II ^ the Sobolev norm 

containing all derivatives up to second order. 

To estabilsh 1^ estimates we use the standard "duality argument." The 

starting point is the following adjoint problem for weV, where the error 

u - u, is the data: 

— ~n 
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(B.2.5) 


T i 

a(v,w) = fv [u-u, JdTt 
- - D _ ^ 


all v£V. 


Suppose for the moment that (B.2.5) can be solved for w_, and 


(B.2.6) II w II 2 < CIIu-i^IIq. 

In this case we put v_ = u_-u^ in (B.2.5) to get 


(B.2.7) = Jlu-UjJ 2 = Hu-u^ll 2 ; 

using orthogonality (i.e* (B.2.1)) we get 
(B.2.8) a(u-H h >w-w h ) = "u-u-Hq 


for any _w n in V, • 


Thus 


(B.2.9) 


II u-u, II llw-w h ll > llu-u, 
u n 


2 

O’ 


We select 


h „h 

Wq E V 


so that 


(B.2.10) 


llw-w h ll < Ch Iw II 2 . 


Thus, with (B.2.6) we get 


(B.2.11) 


D u-u^ II Q < Chll_u-u^H 


Therefore, if linear elements are used 
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ii u-LL n < ch^iiuii 0 . 
n 0 — 2 


The final task is to check (B.2.5) for solvability as well as the a 
priori bound (B.2.6). Rewriting (B.2.5) we get 


(B.2.12) /Lv^LwdTr = Jv^(u-u, )diT all _veV. 

D D 

Suppose vsV and v=0 on T . Then integration by parts gives 


(B.2.13) 

where 


/v^L LwdTT = fv 1 (_u-u, ) d^r 
D D 


all veV, 


* 


L L 


curl curl-grad div 


t. 


Thus defining w_ by 


Lw = _u-u^ 


in D 


w = 0 on T 


it follows that w_ satisfies (B.2.5). Moreover, the a priori bound (B.2.6) 
follows from the theory of second order elliptic equations [8] . 
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Part C 

A Finite Difference Approach 


Cl. Notatlonal Preliminaries; 

Box Variables * 

The notations about to be introduced are most naturally interpreted 
when D can be subdivided into cells {n} each of which is a rectangular 
box. In a later section we shall indicate how more general subdivisions may 
be treated explicitly. 

Following Keller [7] we call v a box variable if it is defined at the 
vertices of cells. For our purpose the importance in employing box variables 
lies in the fact that any average value of a function taken over either a cell 
volume, a face, or an edge of a cell can be approximated in terms of box 
variables by means of the trapezoidal rule. Certain other properties to be 
described then provide a finite difference calculus by means of which 
summation by parts leads to results similar to those established in Section A. 

We employ the notation: v* indicates the centered average and v,^ 

the centered divided difference with respect to x^, i.e., 


V 1 = (v(x i +Ax i /2)+v(x i -Ax i /2))/2 

(C.1.1) 

v ^ = (v(x i +Ax i /2)-v(x i -Ax i /2))/Ax ;j , 
For smooth v, therefore, 


x i +Ax ± /2 


1 = J/ v(x)dx I v Ax 1 + 0(Ax^); 


lX i" Ax d/ 2 


(C. 1.2) 
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in particular, 

(C.1.3) v"^ , = (/9,vdx .dx dx. ) /Air + O(h^) 

, K K 1 j K 


where Att = Ax^Ax^Ax^ and h = min{Ax^}. 


The algebraic identity 


(C. 1.4) 


(vw) 


>i 


i 

= v w 


,1 


i 

w v 




then yields a summation by parts formula while the definition 


v ^ i:j = [' v( x jL +Ax i /2 ,Xj +Ax^ /2 )+v(x i -Ax i /2 -A Xj fl) -v( x 1 -Ax i /2 , Xj +A Xj /2 ) 

-v(x 1 +Ax ± /2,x -Ax : . /2 ) ] * (Ax^x ) 


shows that 
(C.1.5) 

The definitions 
(C. 1.6) 


v . . = v . . 


div^i = l u| k 
i=l ’ 


curl h H = ( u ^ -u ^j ) (i J k) 


provide finite volume approximations to divu and curlu. respectively, i.e. 


div^u = (/divudir) t Att + O(h^) 


curl, u 
1 r 1 


(/curlu.6Tr) v Att + o(h^) . 
7T 


(C • 1 • 7) 
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We propose to examine the finite difference approximations to (1) given, 
in terms of box variables in a cell, by 

.. ijk 

div^u = p 

(C. 1.8a) (ijk) 

curl,u = 
h— — 

The boundary conditions ju*ii = 0 are imposed by 
(C. 1.8b) u^m^ = 0 (ijk) 

where u^ is the trapezoidal approximation to the average value of u on a 
face whose normal is n^. When box variables are understood we shall 

often write ii*ii = 0 to mean the condition expressed by (C.1.8b). 


C2. A Summation Identity 

Define, using box variables, 


(C.2.1) 


t — * 23 . 13 » 12 

L h— = A l— , 1 + A 2— , 2 + A 3— ,3’ 


- _ c 1.23 123> 

Lu = [p >k ) 


where the coefficient matrices are the same as in the definition of Lu_ in 
(A. 1*1). The box-variables scheme (C.1.8) expressing divu » p, curlu = 


may then be written as 



16 


(C.2.2) 

Next, define 

and 

(C.2.3) 


hia-4- 


. - r 23 13 12 ^ 

grad^u = [u ^ ,u }2 ,u >3 J 

(u,v) h = (grad-urgrad-v) . 




The summation-by-parts formula (C . 1 . 4) then leads to 


(C.2.4) 

( L h“) T| 

^L^v) = + 

•u 

where 

is the vector with components 

(C.2.5) 

h 


and 

to 

K 

II 



lc 1c 

Since u is a box-variable u = ii .. ((C.1.5)) and, since 

y ij y J 1 

rp -L 

B = -B. , then ft v = 0* 
k k h— 

Next, multiply (C.2.4) by Att and sum over D using the summation 
analogue of Gauss' theorem to obtain 


(€• 2 . 6 ) 


^(L^^jVAir = I(u,v) h ATT + li^nAix. 


By expressing as given by (C.2.5) in terms of the components of _u 
and v_ (as in (A. 1.7)) the reader may verify that the boundary contribution 
vanishes when u*n — v*n = 0 on T • 


Hence, defining 
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(C.2.7) 


a u(u»v) = l (u,v) 

D h 


Air 


we may state: for any box-variables u_ and v_ satisfying (C.1.8b), i*e», 

ti *ji = v_n = 0 on T, 


(C.2.8) 


a h (u>Z) = 


C3. A Convergence Estimate 

The box-scheme (C.2.2) is an overdetermined system of algebraic equations 

T- 

under the boundary conditions ti*_n = 0 on V . Consider a solution _u_ 11 as 
determined by the least squares problem 

(C.3.1) min I(L h u-f h ) T (L h u-f h )ATT 

with u. *n_ = 0 on T • 

Using (C*2.8) the Euler equations arising from this problem leads to: 
u* provides a least-square solution of if 

(C.3.2) a h (z. h »z) = Il^ L hZ A7r 

for any box-variable jv satisfying v_*n = 0- Write 

(C.3.3) llu h ll h = (a h (u h ,u h )) 1/2 . 

Since only central averages and differences are involved in the 
definition of L^u_ it follows that if u is a solution of Lu_ = which has 
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continuous and bounded mixed third derivatives then 


L^(_u-u^) = O(h^), i.e. 


(C.3.4) llu-u h ll h = 0(h 2 ). 

Define 

(C.3.5) "uHg = [ l l u T (a)u(a)] X/ 2 

ireD aeTT 

where ( cr ) indicates the box-variable approximation to the average value 

of ju_ over a face a of a cell tt. For the continuous problem Friedrichs' 
inequality, when ti*ri = 0, yields 


II ii II Q < y II ii II 


for some constant Y* The same argument which establishes this inequality may 
also be followed using summation and difference operators in place of 
integration and differential operators* The result is 


(C.3.6) 


'a'o 


Y H u II 


Using (C.3.4) we thus obtain 


Hu-u h II q = o(h^), 


i.e. the box-scheme is second order accurate in the II • II q 


norm. 
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C4. A Remark About More General Cells 

As indicated earlier the definitions of div.u, curl.u can be made 

n n — 

independent of any assumption about the shape of a cell n. Properly 
interpreted, so also can the summation by parts formula (C . 1 . 4) as well as 
(C.1.5). A little reflection will convince the reader that the convergence 
proof just given applies as well for irregularly shaped subdomains. 

The followng describes explicit representations for the box-scheme on 
irregular cells. 

Let 0^ denote an oriented face of tt,v = 1,2,..., 6. Applying (C.1.7) to 
a smooth function u_, and employing Gauss' theorem, 

(C.4.1) AiMdiviAi - /divnd^ ~ l f u/dcr. - ii(clJ # A£ 

11 n v a v 

where u.(jl v ) indicates the value of u^ at the centroid of £ v * By 

approximating by the average of its values at the vertex points of 

a an expresison for div.u results when u is a box-variable, 
v n— — 

Similarly, 

curl^u*n^ - /(curl^n^Jdn/An. 

IT 

Apply Stoke 's theorem on a face d£ v and let dsi^ indicate an element of arc 
length in the direction t ^ = do^xii^; the result is 

(C.4.2) curl, u*n t - Y / (u*t .Ids dx, /An = Y(u*t ,)Aa /An 

' hr 1 — i L J v vi ; v i LK vi ; \) 

v a v 

v 

where u/t^ is evaluated at the centroid of the face having the area 

Aa^. By evaluating as the average of its values at the vertices of 

— v (C*4.2) provides an interpretation of curl^u in terms of box-variables . 
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Concluding Remarks 

A comparison of the convergence proofs employed in Parts B and C is of 
interest. The finite-element approach allowed the integral identity (A. 1.9) 
to be used directly but utilized estimates arising from the adjoint problem. 
The finite difference approach, on the other hand, required the development of 
a summation identity corresponding to (A. 1.9); Friedrichs' inequality then 
provided in the required convergence estimate. 

Both proofs are independent of any assumptions about the type of cells 
tt upon which approximations are based. However, the numerical implementation 
of the finite-element approach in such cases may be simpler to employ than the 
box-scheme because of extensive and readily available software for finite- 
element methods. The above remarks suggest that this software could be 
adapted to the finite difference method as well. 

Both methods lead to sparse matrices which may be solved by direct 
algebraic techniques. An iterative scheme for least squares problems due to 
Kaczmarz [6] and Tanabe [11] provides an alternative approach and has been 
employed in [4] to treat the finite difference problem in two dimensions. 
Progress in developing fast iterative methods has been reported to us in 
personal communictions by our colleagues (Grosch and Phillips) and will be 
described elsewhere. 
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